{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interpolation\n",
    "Interpolation is the process of constructing a function that takes on specified values at specified points. In engineering and science, one often has a number of data points, obtained by sampling or experimentation, which represent the values of a function for a limited number of values of the independent variable. It is often required to interpolate (i.e. estimate) the value of that function for an intermediate value of the independent variable.\n",
    "\n",
    "Smile package `smile.interpolation` provides a variety of algorithms on 1d and 2d data. These algorithms implements the interface Interpolation (or Interpolation2D for 2d data), of which the method interpolate takes a value and return an interploated value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import $ivy.`com.github.haifengl::smile-scala:4.0.0`\n",
    "import $ivy.`org.slf4j:slf4j-simple:2.0.16` \n",
    "\n",
    "import java.awt.Color.{BLACK, BLUE, CYAN, DARK_GRAY, GRAY, GREEN, LIGHT_GRAY, MAGENTA, ORANGE, PINK, RED, WHITE, YELLOW}\n",
    "import smile.interpolation._\n",
    "import smile.plot.swing._\n",
    "import smile.plot.show\n",
    "import smile.plot.Render._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Piecewise Linear Interpolation\n",
    "`LinearInterpolation` is quick and easy, but it is not very precise. Another disadvantage is that the interpolant is not differentiable at the control points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "val x = Array(0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0)\n",
    "val y = Array(0.0, 0.8415, 0.9093, 0.1411, -0.7568, -0.9589, -0.2794)\n",
    "\n",
    "val controls = Array.ofDim[Double](x.length, 2)\n",
    "for (i <- 0 until x.length) {\n",
    "  controls(i)(0) = x(i);\n",
    "  controls(i)(1) = y(i);\n",
    "}\n",
    "\n",
    "val linear = new LinearInterpolation(x, y)\n",
    "\n",
    "val data = (0 to 60).map { i =>\n",
    "      val x = i * 0.1\n",
    "      val y = linear.interpolate(x)\n",
    "      Array(x, y)\n",
    "}.toArray\n",
    "\n",
    "val canvas = plot(controls, '*')\n",
    "canvas.add(LinePlot.of(data, Line.Style.SOLID, BLUE, \"Linear\"))\n",
    "show(canvas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`BilinearInterpolation` is an extension of linear interpolation for interpolating functions of two variables on a regular grid. The key idea is to perform linear interpolation first in one direction, and then again in the other direction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cubic Spline Interpolation\n",
    "Spline interpolation uses low-degree polynomials in each of the intervals, and chooses the polynomial pieces such that they fit smoothly together. The resulting function is called a spline. The natural cubic spline is piecewise cubic and twice continuously differentiable. Furthermore, its second derivative is zero at the end points.\n",
    "\n",
    "Like polynomial interpolation, spline interpolation incurs a smaller error than linear interpolation and the interpolant is smoother. However, the interpolant is easier to evaluate than the high-degree polynomials used in polynomial interpolation. It also does not suffer from Runge's phenomenon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "val cubic = new CubicSplineInterpolation1D(x, y)\n",
    "\n",
    "val data = (0 to 60).map { i =>\n",
    "      val x = i * 0.1\n",
    "      val y = cubic.interpolate(x)\n",
    "      Array(x, y)\n",
    "}.toArray\n",
    "\n",
    "canvas.add(LinePlot.of(data, Line.Style.DASH, RED, \"Cubic\"))\n",
    "show(canvas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kriging Interpolation\n",
    "Kriging interpolation can be used for the data points irregularly distributed in space. Kriging belongs to the family of linear least squares estimation algorithms, also known as Gauss-Markov estimation or Gaussian process regression. We implement ordinary kriging for interpolation with power variogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "val kriging = new KrigingInterpolation1D(x, y)\n",
    "val data = Array.ofDim[Double](61, 2)\n",
    "for (i <- 0 to 60) {\n",
    "  data(i)(0) = i * 0.1\n",
    "  data(i)(1) = kriging.interpolate(data(i)(0))\n",
    "}\n",
    "\n",
    "canvas.add(LinePlot.of(data, Line.Style.DOT, CYAN, \"Kriging\"))\n",
    "show(canvas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RBF Interpolation\n",
    "Radial basis function interpolation is a popular method for the data points are irregularly distributed in space. In its basic form, radial basis function interpolation is in the form `y(x) = Σ wi φ(||x-ci||)`, where the approximating function `y(x)` is represented as a sum of N radial basis functions φ, each associated with a different center ci, and weighted by an appropriate coefficient wi. For distance, one usually chooses euclidean distance. The weights wi can be estimated using the matrix methods of linear least squares, because the approximating function is linear in the weights.\n",
    "\n",
    "The points ci often called the centers or collocation points of the RBF interpolant. Note also that the centers ci can be located at arbitrary points in the domain, and do not require a grid. For certain RBF exponential convergence has been shown. Radial basis functions were successfully applied to problems as diverse as computer graphics, neural networks, for the solution of differential equations via collocation methods and many other problems.\n",
    "\n",
    "Other popular choices for φ comprise the Gaussian function and the so called thin plate splines. Thin plate splines result from the solution of a variational problem. The advantage of the thin plate splines is that their conditioning is invariant under scalings. Gaussians, multi-quadrics and inverse multi-quadrics are infinitely smooth and and involve a scale or shape parameter, r0 > 0. Decreasing r0 tends to flatten the basis function. For a given function, the quality of approximation may strongly depend on this parameter. In particular, increasing r0 has the effect of better conditioning (the separation distance of the scaled points increases).\n",
    "\n",
    "A variant on RBF interpolation is normalized radial basis function (NRBF) interpolation, in which we require the sum of the basis functions to be unity. NRBF arises more naturally from a Bayesian statistical perspective. However, there is no evidence that either the NRBF method is consistently superior to the RBF method, or vice versa."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "val rbf = new RBFInterpolation1D(x, y, new smile.math.rbf.GaussianRadialBasis())\n",
    "val data = Array.ofDim[Double](61, 2)\n",
    "for (i <- 0 to 60) {\n",
    "  data(i)(0) = i * 0.1\n",
    "  data(i)(1) = rbf.interpolate(data(i)(0))\n",
    "}\n",
    "\n",
    "canvas.add(LinePlot.of(data, Line.Style.DOT_DASH, GREEN, \"RBF\"))\n",
    "show(canvas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Shepard Interpolation\n",
    "Shepard interpolation is a special case of normalized radial basis function interpolation if the function φ(r) goes to infinity as `r → 0`, and is finite for `r > 0`. In this case, the weights wi are just equal to the respective function values yi. So we need not solve linear equations and thus it works for very large N. An example of such φ is `φ(r) = r^-p` with (typically) `1 < p ≤ 3`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "val shepard = new ShepardInterpolation1D(x, y, 3)\n",
    "val data = Array.ofDim[Double](61, 2)\n",
    "for (i <- 0 to 60) {\n",
    "  data(i)(0) = i * 0.1\n",
    "  data(i)(1) = shepard.interpolate(data(i)(0))\n",
    "}\n",
    "\n",
    "canvas.add(LinePlot.of(data, Line.Style.LONG_DASH, PINK, \"Shepard\"))\n",
    "show(canvas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Shepard interpolation is rarely as accurate as the well-tuned application of other radial basis functions. However, it is simple, fast, and often sufficient for quick and dirty applications."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2D Interpolation\n",
    "\n",
    "For 2d grid, we provide two implementations: `CubicSplineInterpolation2D` and `BicubicInterpolation`. `CubicSplineInterpolation2D` is similar to one-dimensional splines as it guarantees the continuity of the first and second function derivatives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "val x1 = Array(1950.0, 1960, 1970, 1980, 1990)\n",
    "val x2 = Array(10.0, 20, 30)\n",
    "val y = Array(\n",
    "            Array(150.697, 199.592, 187.625),\n",
    "            Array(179.323, 195.072, 250.287),\n",
    "            Array(203.212, 179.092, 322.767),\n",
    "            Array(226.505, 153.706, 426.730),\n",
    "            Array(249.633, 120.281, 598.243)\n",
    "        )\n",
    "\n",
    "val canvas = heatmap(y, Palette.jet(256))\n",
    "canvas.setTitle(\"Original\")\n",
    "show(canvas)\n",
    "\n",
    "val cubic = new CubicSplineInterpolation2D(x1, x2, y)\n",
    "\n",
    "val data = Array.ofDim[Double](101, 101)\n",
    "for (i <- 0 to 100; j <- 0 to 100)\n",
    "    data(i)(j) = cubic.interpolate(1950 + i*0.4, 10 + j*0.2)\n",
    "\n",
    "val cubicPlot = heatmap(data, Palette.jet(256))\n",
    "cubicPlot.setTitle(\"Cubic\")\n",
    "show(cubicPlot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In contrast, `BicubicInterpolation` guarantees continuity of only gradient and cross-derivative. Second derivatives could be discontinuous."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "val bicubic = new BicubicInterpolation(x1, x2, y)\n",
    "\n",
    "for (i <- 0 to 100; j <- 0 to 100)\n",
    "    data(i)(j) = bicubic.interpolate(1950 + i*0.4, 10 + j*0.2)\n",
    "\n",
    "val bicubicPlot = heatmap(data, Palette.jet(256))\n",
    "bicubicPlot.setTitle(\"Bicubic\")\n",
    "show(bicubicPlot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In image processing, bicubic interpolation is often chosen over bilinear interpolation or nearest neighbor in image resampling, when speed is not an issue. Images resampled with bicubic interpolation are smoother and have fewer interpolation artifacts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Laplace Interpolation\n",
    "Laplace interpolation can restore missing or unmeasured values on a 2-dimensional evenly spaced regular grid. In some sense, Laplace interpolation produces the smoothest possible interpolant, which are obtained by solving a very sparse linear equations with biconjugate gradient method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "val zz = Array.ofDim[Double](101, 101)\n",
    "val ww = Array.ofDim[Double](101, 101)\n",
    "for (i <- 0 to 100; j <- 0 to 100) {\n",
    "  zz(i)(j) = if (java.lang.Math.random() < 0.2) Double.NaN else data(i)(j)\n",
    "  ww(i)(j) = zz(i)(j)\n",
    "}\n",
    "\n",
    "val canvas = heatmap(ww, Palette.jet(256))\n",
    "canvas.setTitle(\"Missing Values\")\n",
    "show(canvas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LaplaceInterpolation.interpolate(zz)\n",
    "val canvas = heatmap(zz, Palette.jet(256))\n",
    "canvas.setTitle(\"Laplace\")\n",
    "show(canvas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Scala (2.13)",
   "language": "scala",
   "name": "scala213"
  },
  "language_info": {
   "codemirror_mode": "text/x-scala",
   "file_extension": ".scala",
   "mimetype": "text/x-scala",
   "name": "scala",
   "nbconvert_exporter": "script",
   "version": "2.13.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
